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Abstract 

We use a generic and general numerical method to obtain solutions for 
the flow of generalized Newtonian fluids through circular pipes and plane 
slits. The method, which is simple and robust can produce highly accurate 
solutions which virtually match any analytical solutions. The method is based 
on employing the stress, as a function of the pipe radius or slit thickness 
dimension, combined with the rate of strain function as represented by the 
fluid rheological constitutive relation that correlates the rate of strain to stress. 
Nine types of generalized Newtonian fluids are tested in this investigation and 
the solutions obtained from the generic method are compared to the analytical 
solutions which are obtained from the Weissenberg-Rabinowitsch-Mooney- 
Schofield method. Very good agreement was obtained in all the investigated 
cases. All the required quantities of the flow which include local viscosity, 
rate of strain, flow velocity profile and volumetric flow rate, as well as shear 
stress, can be obtained from the generic method. This is an advantage as 
compared to some traditional methods which only produce some of these 
quantities. The method is also superior to the numerical meshing techniques 
which may be used for resolving the flow in these systems. The method is 
particularly useful when analytical solutions are not available or when the 
available analytical solutions do not yield all the flow parameters. 

Keywords: fluid dynamics; rheology; pipe flow; slit flow; Newtonian; power 
law; Ellis; Ree-Eyring; Carreau; Cross; Bingham; Herschel-Bulkley; Casson. 
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1 Introduction 


There are many applications for the flow of generalized Newtonian fluids in circular 
pipes and plane slits. These include biological flow in living organisms and flnid 
shipping in transport and process industries. Hence various analytical and numerical 
methods have been developed and used to obtain the flow parameters which include 
the stress function, rate of shear strain, flow velocity profile, local viscosity and 
volumetric flow rate as functions of the condnit geometry, pressure drop and flnid 
rheology [1, 2]. 

A very simple numerical method, which is very natural to use, is to employ the 
stress function, which is defined in terms of the condnit geometry and pressnre 
drop, directly. The method is based on a simple combination of the easily obtained 
stress fnnction with the flnid rheology to obtain the rate of strain as a fnnction of 
the velocity-varying spatial dimension of the condnit, i.e. the radins for the pipe 
flow and the thickness dimension for the slit flow. The flow velocity profile and the 
volumetric flow rate can then be obtained from simple numerical integrations. The 
local viscosity can also be obtained from the flnid rheological relation as soon as 
the stress and rate of strain are obtained. 

The method is very generic, as well as being more convenient to implement 
and nse, and hence the solntions obtained from it can be more reliable than the 
solntions obtained from more sophisticated methods. It is also very general since 
it can be applied whenever the rate of strain can be expressed as an explicit or 
implicit function of the shear stress; a relation that usually can be obtained from the 
flnid rheological eqnation. Hence, the method can be used to solve almost all the 
generalized Newtonian flow problems in pipes and slits since the stress function for 
these systems can be easily obtained and the appropriate rheological equations can 
be readily expressed in the desired form. The method is particnlarly nsefnl when 
the flow parameters, or some of which, cannot be obtained analytically. Hence, it 
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can provide a good alternative to the classical methods for solving these problems 
which are traditionally solved by employing rather complicated methods such as 
numerical discretization techniques like finite element and finite difference. 

In section § 2 we give a general description of the method and its theoretical 
justification as well as the type of flow systems to which the method applies and 
the restrictions and assumptions that should be observed in its application. In 
section § 3 we discuss practical issues about the implementation of the method. 
We also present a sample of the volumetric flow rate solutions that were obtained 
from the method with comparison to similar results obtained from the Weissenberg- 
Rabinowitsch-Mooney-Schofield (WRMS) method. The paper is ended in section § 4 
with general discussions and summary of the main objectives and accomplishments 
of this investigation. 
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2 Method 


We assume an incompressible, laminar, isothermal, steady, pressure-driven, fully- 
developed flow of purely-viscous, time-independent fluids that are properly charac¬ 
terized by the following generalized Newtonian fluid model 


r = ( 1 ) 

where the shear viscosity, fi, and shear stress, r, depend only on the contemporary 
rate of strain, 7 , with no memory for the fluid of its deformation history. The effects 
of external body forces, like gravitational attraction or electromagnetic interaction, 
as well as the edge effects at the entry and exit zones of the conduit are assumed 
insignificant. Dependencies on physical factors like temperature, which are not 
related to deformation, are also ignored assuming fixed conditions or negligible 
contribution from these factors. The flow is also assumed to be in shear mode with 
no significant extensional contributions. 

As for the boundary conditions, a no-slip at the conduit wall is assumed and 
hence a zero velocity condition at the fluid-solid interface is maintained. For the 
investigated types of rheology and flow system, the flow velocity profile has a 
stationary derivative point at the symmetry center line of the pipe and symmetry 
center plane of the slit which means zero stress and rate of strain at these loci. For 
viscoplastic fluids, these stationary zones extend to include all the points at the 
forefront of the flow profile whose stress falls below the fluid yield stress. 

Concerning the type of conduit, we use circular pipe geometry with length L and 
radius r, and plane slit geometry with length L, width W and thickness 2B where 
the latter dimension is the smallest of the three. A pressure drop Ap is imposed 
along the conduit length dimension which defines the flow direction. The pipe is 
assumed straight with a cross sectional area that is uniform in shape and size while 
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the slit is assumed straight, long and thin with a uniform cross section. For both 
types of condnit, rigid mechanical properties of the condnit wall are assnmed and 
hence the conduit wall is not deformable under the considered range of pressure 
drop. It is also assumed that the slit is positioned symmetrically in its thickness 
dimension, z, with respect to the plane z = 0. 

The generic method is based on the proportionality between the stress and the 
spatial coordinate in the flow profile dimension, s, which stands for the radius r 
in the case of pipe and for the thickness dimension 2 ; in the case of slit. Since the 
stress r as a function of s is known from the above mentioned proportionality, then 
the rate of strain, 7 , can be easily obtained from the flnid rheological constitntive 
relation which correlates the rate of shear strain to the shear stress as long as it 
can be put in the form 7 = 7 (r) where the dependency of 7 on r can be explicit 
or implicit. If 7 is an explicit function of r, as it is the case for example for Ellis 
fluids (refer to Table 1), then 7 can be obtained directly by a simple substitution in 
the rheological relation. If, on the other hand, 7 is an implicit fnnction of r, as it 
is the case for example for Cross flnids (refer to Table 1), then 7 can be obtained 
numerically using a simple numerical solver based for instance on a bisection method. 
In both cases, the obtained rate of strain as a function of s can be used to obtain the 
velocity profile and subsequently the volumetric flow rate by consecutive numerical 
integrations. 

To be more specific, the stress in the investigated flow systems of circular pipes 
and plane slits is a function of the conduit geometry and pressure drop. The shear 
stress as a function of the velocity-varying dimension s can be obtained by several 
methods which can be found in many fluid mechanics and rheology textbooks, 
e.g. [1-3], and hence we are not going to give a formal derivation; instead we give 
the final results with a sketch of how it can be obtained. From a simple force 
balance argument based on applying the Newton second law of mechanics to the 
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non-accelerating steady flow, where the normal force exerted on the condnit cross 
section is balanced by the shear stress force, it can be established that the stress 
is proportional to s as summarized in the following equations for pipe and slit 
respectively 



and 



( 2 ) 


where tr and tr are, respectively, the pipe and slit wall shear stress as given by 


RAp BAp 

and “ 

where R is the pipe radius, B is the slit half thickness, and L is the conduit length 
across which a pressure drop Ap is exerted. 

Since the spatial dependence of the stress function is known, the rate of strain 
as a function of s can be obtained from the rheological relation expressed in the 
form 7(r(s)), as given for a sample of generalized Newtonian fluids in Table 1 , 
and hence 7(5) is obtained. The obtained strain rate function is then integrated 
numerically with respect to the spatial coordinate of the flow profile, s, to obtain the 
flow velocity profile where the no-slip at wall condition provides an initial value for 
the flow velocity, n = 0, which is then incremented in moving from the wall to the 
center during the integration process. The numerically obtained flow velocity is then 
integrated numerically with respect to the conduit cross sectional area perpendicular 
to the flow direction to obtain the volumetric flow rate. For viscoplastic fluids, the 
zero stress condition at the pipe center line and slit center plane is extended to 
include all the regions at the forefront of the flow profile whose shear stress falls 
below the fluid yield stress. In Table 1 the rheological constitutive relations for nine 
fluid models employed in this study are presented. For Carreau and Cross models, 
7 is given as an implicit function of r and hence a simple numerical solver like 
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bisection is required to obtain 7 as a function of r(s) and hence as a function of s. 


Table 1 : The rate of shear strain, 7, as a function of shear stress, r, for the nine 
fluid models used in this investigation [ 1 , 2 , 4 - 6 ]. For Carreau and Cross models, 7 
is given as an implicit function of r. The meaning of the symbols can be obtained 
from Nomenclature S 5 . 


Model 


Rate of Shear Strain 


Newtonian 7 

Power Law 7 

Ellis 7 

Ree-Eyring 7 

Carreau 7 

Cross 7 

Bingham 7 

Herschel-Bulkley 7 

Casson 7 


1^0 




1 + T 


OL —1 


= ^ sinh — 

fir y Tc 


hi + (ho - hi) (1 + 

,, I Mo-Mi 
U'l “T 


= r 


1 +A’"7” 


= r 


_ t-tq 

c 


_ nj T—TQ 

K 


7 















3 Implementation, Results and Assessment 


The generic method, as described in the last section, was implemented in a computer 
code using standard numerical integration and bisection solver techniques. The 
method was then employed to obtain solutions for the flow of nine types of generalized 
Newtonian fluid through circular pipes and plane slits. The nine types of fluid 
are: Newtonian, power law, Ellis, Ree-Eyring, Carreau, Cross, Bingham, Herschel- 
Bulkley and Casson. The numerical results of the volumetric flow rate which are 
obtained from the generic method using wide ranges of fluid and conduit parameters 
were thoroughly compared to the results obtained from the WRMS analytical 
method. 

A representative sample of the results obtained from the two methods are 
presented in Eigures 1 and 2 where the fluid and conduit parameters of these 
examples are given in Table 2. As seen in these examples, the solutions obtained 
from the generic method agree very well with the WRMS analytical solutions. 
The minor departure between the two methods is due mainly to the nature of the 
generic method as it heavily relies on numerical techniques, i.e. bisection solvers 
and successive numerical integrations, which can accumulate numerical errors. 

The big advantage of using the generic method in the flow problems through 
circular pipes and plane slits is that it is very general as it can be applied to 
almost any generalized Newtonian fluid model regardless of the complexity of its 
constitutive relation as long as the rheological relation can be casted as an explicit 
or implicit function of the form 7 (t(s)). Moreover, as the method only involves 
very simple and reliable numerical techniques, like bisection solvers and numerical 
integration techniques, high accuracy of the solution can be achieved. R also has 
no convergence difficulties under normal conditions, as it does not require matrix 
solvers and similar sophisticated computational techniques, and hence it can be 
very reliable. Eurthermore, the method is very simple and hence it is easy and 
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convenient to implement and nse; which is an advantage on its own pins being a 
farther contribnting factor to the reliability of the solntions. Hence, the generic 
method can prodnce virtnally exact solntions that match in their accnracy any 
actnal or potential analytical solntions even for the flnids with complex rheology. 
The generic method can therefore offer a better alternative to the commonly used 
numerical methods which employ more sophisticated computational techniques like 
finite difference. 

Although the generic method is presented here as a numerical technique, it 
can also be applied analytically where the stress function, r(s), is substituted in 
the rheological relation to obtain the rate of strain, 7(5), which can be integrated 
successively to obtain the flow velocity profile and volumetric flow rate. This can 
provide alternative analytical forms that can be useful in some cases for verification 
or other purposes. 


Table 2 : Fluid and conduit parameters for the examples of Figures 1 and 2 . The 
‘R’ column applies to pipes and the ‘B’ column applies to slits; the other columns 
are common to both. SI units apply to all dimensional quantities as given in 
Nomenclature § 5 , and ‘HB’ stands for Herschel-Bulkley. 


Model 

Flnid Properties 

R 

B 

L 

Newtonian 

Po = 0.013 

0.05 

0.007 

0.45 

Power Law 

k = 0.023, n = 0.63 

0.015 

0.008 

0.55 

Ellis 

Pe = 0.026, Th = 8, a = 1.6 

0.005 

0.002 

0.22 

Ree-Eyring 

CO 

T—1 

0 

II 

0.03 

0.004 

1.65 

Carrean 

po = 0.37, Pi = 0.0079, A = 0.65, n = 0.62 

0.043 

0.0086 

0.95 

Cross 

po = 0.42, Pi = 0.0093, A = 0.77, m = 0.58 

0.032 

0.01 

1.26 

Bingham 

C = 0.033, TO = 5.7 

0.01 

0.01 

0.25 

HB 

C = 0.042, To = 7.9, n = 1.34 

0.005 

0.002 

0.13 

Casson 

K = 0.71, To = 3.2 

0.08 

0.011 

1.27 
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(a) Newtonian (b) Power Law (c) Ellis 




Figure 1: Comparing the WRMS analytical solutions (solid line) with the numerical 
solutions from the generic method (circles) of Q in m^.s“^ (vertical axis) versus Ap 
in Pa (horizontal axis) for the flow of the nine fluid models in pipes. The pipe and 
fluid parameters are given in Table 2. 
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(a) Newtonian 


(b) Power Law 


(d) Ree-Eyring 


(e) Carreau 


(c) Ellis 


(f) Cross 


(g) Bingham 


(h) Herschel-Bulkley 


(i) Casson 




Figure 2: Comparing the WRMS analytical solutions (solid line) with the numerical 
solutions from the generic method (circles) of Q in m^.s“^ (vertical axis) versus Ap 
in Pa (horizontal axis) for the flow of the nine flnid models in slits. The slit and 
flnid parameters are given in Table 2. In all these examples W = 1.0 m. 
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4 Conclusions 


In this study we presented a generic and general numerical method for finding the 
flow solutions of generalized Newtonian fluids in one dimensional flow problems 
that can be applied easily to circular pipes and plane slits. The method can be 
used to obtain all the required flow parameters which include shear stress, local 
viscosity, shear rate, flow velocity profile and volumetric flow rate. 

Thorough comparisons were made between the results of the generic method and 
the results of the analytical solutions obtained from the Weissenberg-Rabinowitsch- 
Mooney-Schofield method. In all cases, the two methods produced virtually identical 
results considering the numerical errors introduced by the heavy use of numerical 
methods, like numerical integration and bisection solvers, in the generic method. 

The generic method enjoys several advantages over the competing numerical 
methods like finite element and finite difference. These advantages include ease and 
convenience to implement and use, generality, reliability and accuracy. The method 
may also be useful to apply analytically in some cases. 

The generic method can be particularly useful when no analytical solutions 
can be obtained, or the analytical solutions can provide only some of the flow 
parameters, e.g. volumetric flow rate, but not others, e.g. flow velocity profile. 
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Nomenclature 


B slit half thickness (m) 

C viscosity coefficient in Herschel-Bulkley model (Pa.s"^) 

C viscosity coefficient in Bingham model (Pa.s) 
k viscosity coefficient in power law model (Pa.s"^) 

K viscosity coefficient in Casson model (Pa.s) 

L conduit length (m) 

m indicial parameter in Cross model 

n flow behavior index in power law, Carreau and Herschel-Bulkley models 

Ap pressure drop across the conduit length (Pa) 

Q volumetric flow rate (m^.s“^) 
r radius (m) 

R pipe radius (m) 

s spatial coordinate that represents r for pipe and for slit (m) 

V fluid velocity in the flow direction (m.s“^) 

W slit width (m) 

spatial coordinate of slit thickness dimension (m) 

a indicial parameter in Ellis model 
7 rate of shear strain (s“^) 

A characteristic time constant in Carreau and Cross models (s) 
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jj, fluid shear viscosity (Pa.s) 

/io zero-shear viscosity in Carreau and Cross models (Pa.s) 

yU-e low-shear viscosity in Ellis model (Pa.s) 

fii infinite-shear viscosity in Carreau and Cross models (Pa.s) 

/io Newtonian viscosity (Pa.s) 

fir characteristic viscosity in Ree-Eyring model (Pa.s) 
r shear stress (Pa) 

To yield stress in Bingham, Herschel-Bulkley and Casson models (Pa) 
tb shear stress at slit wall (Pa) 

Tc characteristic shear stress in Ree-Eyring model (Pa) 

Th shear stress when viscosity equals ^ Ellis model (Pa) 
tr shear stress at pipe wall (Pa) 
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